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Abstract. We have adapted the maximally-localized Wannier function approach of 
[Souza I, Marzari N and Vanderbilt D 2002 Phys. Rev. B 65 035109] to the density 
functional theory based Siesta code [Soler J M et al. 2002 J. Phys.: Cond. Mat. 14 
2745] and applied it to the study of Co substitutional impurities in bulk copper as well 
as to the Cu (111) surface. In the Co impurity case, we have reduced the problem 
to the Co d-electrons and the Cu sp-band, permitting us to obtain an Anderson-like 
Hamiltonian from well defined density functional parameters in a fully orthonormal 
basis set. In order to test the quality of the Wannier approach to surfaces, we have 
studied the electronic structure of the Cu (111) surface by again transforming the 
density functional problem into the Wannier representation. An excellent description 
of the Shockley surface state is attained, permitting us to be confident in the application 
of this method to future studies of magnetic adsorbates in the presence of an extended 
surface state. 
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1. Introduction 

The advent of density functional theory (DFT) in electronic structure calculations has 
revolutionized the fields of quantum chemistry and, more generally, of condensed matter 
physics pp. However, due to necessary approximations of the unknown functionals 
(typically, the local density approximation or LDA, and semilocal approximations as 
the generalized gradient approximation or GGA), important limitations prevail in the 
description of many systems ranging from insulating materials to impurities in a metal 
host. The general approach to these problems is to go beyond DFT by proposing a 
simplified Hamiltonian that can be solved, at the expense of losing the parameter-free 
advantage of DFT. Recently, the situation is changing and new methods either stemming 
from DFT or making use of DFT for initial input are emerging [2j [3]. One such example 
is the LDA+U approach [U [5] where the exchange-correlation potential for the local 
electron gas is complemented by missing strong localized correlations. The success of 
this approach has been considerable in explaining the opening of a band gap in systems 
with localized electronic states. Yet, the choice of the Coulomb intra-atomic energy U 
is somewhat arbitrary, depending on the choice of basis sets or other descriptors of the 
treated system [6]. 

Nakamura and co-workers [7J have developed a constrained LDA approach based 
on maximally localized Wannier functions (MLWF) |8j. The MLWF replace the linear 
muffin-tin orbitals (LMTO) that were initially used as a natural basis set to define 
U [HEIDI]]. As Nakamura et al. emphasize [7J, the properties of MLWF are particularly 
appealing for reducing the complicated DFT problem to a simplified Hamiltonian, where 
especial physics can be explored such as the localized electron correlations mentioned 
above. Wannier functions have been thoroughly explored in the work by Marzari 
and Vanderbilt [UJ, and an algorithmic approach to obtaining them from a DFT 
calculation is available. More recently, Souza, Marzari and Vanderbilt [8] have extended 
the approach to disentangle electronic bands, creating a compact local basis set that 
accurately reproduces the DFT electronic structure in a given energy window. This 
approach is independent of the actual implementation of the DFT calculation, yielding 
a natural way of describing an extended basis set calculation in terms of localized 
functions. 

In the present work, we have interfaced the approach by Souza et al. [8] as 
implemented in Wannier90 [12] to the Siesta code [13]. Contrary to the case of 
Nakamura et al. [7J, Siesta is an atomic basis set code, and it would seem natural 
to use the atomic orbitals to define U and associated model Hamiltonians. However, 
besides their optimized spread, MLWF have two important features: (i) they are a 
naturally orthogonal basis set, rendering tight-binding like approaches easy to use (ii) 
the extraordinary accuracy of MLWF in a defined energy window permits to disentangle 
electronic bands [8] and to have a simple tight-binding approach with DFT accuracy 
within the chosen energy window. Hence, our present implementation of MLWF permits 
us to translate the complex atomic orbital method into a simple orthogonal tight-binding 
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approach that can be easily cast into an Anderson-Hamiltonian model [H] or a Hubbard 
one [15] . 

The manuscript is organized in four sections. The methods section deals with 
a small overview of the work of reference [8], and the technical details of the actual 
numerical implementation. The approach is then applied to two different systems: (i) 
a substitutional Co impurity in bulk Cu and (ii) the Cu (111) surface. Both systems 
are of uttermost interest. The case of Co in Cu is a classical example where strong 
correlations are in play leading to very large Kondo temperatures [161 E]> hence it is 
interesting to understand the magnetic properties leading to the Kondo correlations in 
this system [18J. Besides, Co in Cu is one of the model giant magnetoresistance (GMR) 
systems [191 EQ1 [21] . The concentration of Co in Cu is an important parameter in its 
magnetic properties [21]. Indeed, when the Co to Cu atomic ratio is above 1 to 4, the 
system becomes a ferromagnet. Below this density, Fan et al. (21] experimentally find 
that the Co atoms are distant enough to show paramagnetism in agreement with the 
paradigmatic analysis by Goodenough [22]. We will concentrate in this paramagnetic 
phase, with two different densities, for 2x2x2 and 4x4x4 cubic supercells. 

Our calculations are a first step in the study of the complex electronic structure 
of the Co-impurity problem that presents a high (~ 500 K) Kondo temperature |TT] . 
Hence, the evolution of the electronic structure as the Co concentration is reduced, 
permits us to study the effect of Co concentration and the use of MWLF gives us access 
to an Anderson-like Hamiltonian. 

The study of the Cu(lll) surface is also very interesting given the existence of the 
LL' gap (f in the surface Brillouin zone) and the associated Shockley state. Despite the 
locality of the MLWF, the calculations succeed in accounting for the surface Shockley 
state and in yielding the correct electronic structure about the Fermi energy. 

This work shows that MLWF are accurate enough to study the magnetic properties 
of Co on Cu(lll) and are a first step towards the study of the electronic structure [23] 
on surface systems as well as the Kondo effect of magnetic adsorbates that has recently 
received much experimental [211 EH] and theoretical [26] attention. 

2. Method 

2.1. Wannier functions and Wannier90 

Wannier functions, w(r — R), are formally identical to Fourier coefficients of Bloch 
waves, V'k(r), given by 



Here R denote Bravais vectors and Q* is the volume of the Brillouine zone (BZ). 
This definition suffers from the indeterminacy of the phases of Bloch functions, ^k(r). 
Furthermore, if a group of bands, {n} is considered, with corresponding Bloch functions, 
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^kn( r ); it is ambiguous to talk about an individual band and the above definition can 
be generalized to 

w m (r-K) = ^ I ^f/ m „(k)Vkn(r)e- ik - R d 3 k, (2) 

where U mn (k) is an arbitrary unitary matrix. The unitarity guarantees that the Wannier 
functions will be orthogonal. The arbitrariness of U mn (k) allows for tuning the phases 
of Bloch functions in the integral as well as the admixture of functions pertaining to 
different bands. Thus, there is a whole class of Wannier functions for the given band 
structure. Marzari and Vanderbilt devised a variational scheme that determines the 
Wannier functions with minimum total spread, Q defined as 

" = £ ^ - < r >™] > ( 3 ) 

m 

where we have used the notation ( • ) m = (w m \ ■ \w m ). For a given set of Bloch functions, 
the total spread, Q, is a functional of the unitary matrices U mn (\t). The Wannier 
functions so obtained are called maximally- localized Wannier functions (MLWF). 

The locality and the orthogonality of Wannier functions permit us to have a 
straightforward compact tight-binding representation of the Hamiltonian in systems 
where the group of bands of interest is separated from the rest of the electronic structure 
by a band gap. In metals, the absence of band gaps renders the separation of states more 
complicated. Souza et al. have devised a disentanglement procedure [H] by focusing on a 
certain energy window, hereinafter called outer energy window, and by selecting certain 
bands in it. Hence, the disentanglement procedure necessarily reduces the number of 
states inside the outer energy window, while exactly reproducing the electronic bands in 
a certain energy interval called inner energy window. The selection of bands proceeds 
via trial orbitals g n (r) that allow to define the character of the Wannier subspace of 
interest. 

The numerical calculations of MLWF have been done with the Wannier90 [12J 
code. This package can be used as a post-processing tool with most first-principles codes. 
In practical simulations, the Bloch states, \ipim), are computed in a mesh of uniformly 
spaced k-points within the first-Brillouin zone. The basis set used to represent the 
wave functions changes from one electronic structure code to another one; however, the 
MLWF algorithm requires an input that is essentially basis- independent. In particular, 
the main ingredients are (i) the overlap matrix between the cell-periodic parts of the 
wavefunctions at neighboring k-points 

M mn (k, b) = (u km \u k+hn ) = (^km|e" lb ' r |Vk+bn), (4) 

(ii) the Bloch energies £kn on the regular grid of k-points and (iii) the coefficients of the 
above trial orbitals, g n {v), in the Bloch basis 



■^mn(k) — \lpkm\9n) ■ 



(5) 
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The trial functions we employ are of the form g n (r) = R n (r)Qi mr ((p,9), r = r,4>,6, 
a product of nodeless hydrogenic-like radial part R n (r) and a real spherical harmonic 
®im r ((p,0) with angular momentum I and projection m r . 

From this, the Wannier90 code computes the final unitary transformation 
matrices C/ TOn (k). With the aid of the matrices U mn (k) the Hamiltonian can be expressed 
in the Wannier basis and diagonalized at any k-point. These Wannier-based energy 
bands are further compared to the initial ab-initio bands in order to test the quality of 
the newly obtained Wannier basis set. 

It is very useful to work with density of states projected onto a certain Wannier 
function of spin cr, w ma . This projected density of states (PDOS)[27| is the spectral 
function, p ma (u), given by 

Pma{u) = ^2\(w ma \?p kna )\ 2 5(U - 6 kna ) 
nk 

= ^2 \ u mn(ka)\ 2 5(u - e kn<7 ). (6) 

nk 

This Wannier PDOS or spectral function shows how the Wannier character is distributed 
in energy over the electronic band structure 0. 

From ([6]) it is straightforward to define the Wannier function occupation at zero 
temperature as 

n ma = / p ma (u)du, (7) 

J — oo 

where /i is the Fermi level. 



2.2. Pseudo-atomic orbital DFT calculations 

Ab-initio DFT calculation presented in this work are based on strictly localized [28J 
numerical pseudo-atomic orbitals (PAO) [29] that are solutions to the atomic Kohn and 
Sham equation with norm- conserving pseudopotentials. In particular, the calculations 
are done using the Siesta [13J package. 

The matrix elements between basis functions are calculated by real-space 
integration, and the Hamiltonian eigenstates are labeled using the Bloch theorem, 
because periodic boundary conditions are imposed. The Bloch functions can be 
expanded into the PAO basis as follows 

VVk(r) = c M n(k)e* k - (l > +R VM(r - r„ - R) (8) 

where we assume that the unit cell contains the centers r M of basis functions </v( r — 
which are then repeated periodically to every other cell by the Bravais lattice vector R. 

| The PDOS can be used to analyze any electronic structure by choosing the analyzing functions, for 
example, molecular orbitals were used to describe the electronic structure of benzene on Cu(100) in 

E3- 
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Finally, the complex numbers c Mn (k) are expansion coefficients. The Brillouin zone is 
sampled uniformly. For the exchange and correlation potential entering the Kohn and 
Sham equation, we use the PBE generalized gradient approximation [30]. The chosen 
atomic basis set is an optimized double-^ plus polarization for the valence states of Co 
and Cu, amounting to 15 basis functions per atom. All the parameters that define the 
shape and the range of the basis functions were obtained by a variational optimization 
of the enthalpy (energy plus a penalty for orbital volume increase) with a pressure 
of P = 0.1 GPa, following the recipe given in reference [31]. The substitutional Co 
impurity basis set was optimized in the Cu host, and the Cu basis set corresponds to 
the optimal one for bulk Cu. The question of optimal basis sets for surfaces is more 
intricate and has been discussed in [32]. We use the basis set obtained in reference [32J, 
which is based on the same enthalpy minimization with an especial focus on the vacuum 
extension of the surface state density. 



2. 3. Implementation of Wannier functions 

The implementation of maximally localized Wannier functions in Siesta consists in 
evaluation of (j3J) and (jSJ). Expanding the Bloch functions according to (jHJ), we obtain 

M mn (k,b) =^^ C ; m (k)c, n (k + b)e ik -( R - r ^M^(R,b) (9) 

iiv R 



along with 



where 



and 



A mn (k) = ^c; m (k) ^e- k -( R+r ^^ n (R), (10) 



R 



M^(R, b) = / <p*(r + R - iv + r,)e- Jb - r ^(r)d d r (11) 



^ n (R) = J ^(r - r„ - R)s„(r)d 3 r. (12) 

This reduces the computation of (j3J) and (jSJ) to the calculation of a few matrix 
elements of localized functions, equations ( flTl) and (1T2"]) . The first integral is computed 
on the real space grid while the second integral makes use of the analytic angular 
dependence of the integrand in the same way as is done for the calculation of overlap 
matrices in reference |13j. This is an important difference with the implementation of 
[33] where an expansion on powers of the integrand is performed. Another difference of 
our implementation is that we write an interface for Wannier90, and hence the trial 
functions are the ones for Wannier90, while in reference [33] the original basis set is 
used. 

Finally, we note that Brillouin zone sampling used to obtain the self-consistent 
electronic density in Siesta is essentially independent from the sampling used for the 
input data of Wannier90, i.e. Equations ( Hf5l) and e n k- 
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2.4- Ab-initio calculation of a cobalt impurity in bulk copper 

The first system, the cobalt impurity in a bulk FCC copper host matrix, was simulated 
by a supercell where one host atom was replaced by a cobalt one. The lattice parameter 
was fixed to the theoretical value we found for pure copper (a/ at = 3.690A). We use two 
different unit cells to represent the system: (i) an eight-atom (2x2x2) cell [§| , where 
there is a 1/7 ratio of Co impurities, and (ii) a 64- atom one (4 x 4 x 4) [jj where the cobalt 
concentration drops to 1/63. Besides its lower computational demand, the small cell 
permits us to represent the calculated bands and follow the Wannier disentanglement 
in a simpler way avoiding the larger band folding of the (4x4x4) supercell. According 
to the experimental data by Fan et al. [2T] only above a concentration ratio of 1/4 does 
Co in Cu show ferromagnetic ordering; below this concentration, the system becomes 
paramagnetic. Our calculations correspond to impurity densities of the paramagnetic 
phase. 

Interestingly, the measured saturation magnetic moment per atom is 1.5 fis for 
pure Co and it drops to 0.4/is/atom for a 1/9 concentration [21J. Our calculations 
yield 0.15yUs/atom for the 1/7 concentration, indicating a smaller spin-polarization of 
our DFT calculations. It is difficult to know the actual experimental error bar, but 
our calculations yield the right order of magnitude as well as the trend with decreasing 
Co concentration. Indeed, the supercell magnetization for the 2 x 2 x 2-case (1/7 
concentration) is 1.03 /i# (0.15yUs/atom) while the supercell magnetization for the 
4x4x4 is 1.57 fiB (0.025/iB/atom), indicating a saturation of the cell magnetization and 
thus a decrease of magnetization per atom, roughly linear with the number of electrons 
following the Slater-Pauling curve [21]. 

For the obtention of the MLWF, 55 bands in the 2x2x2 supercell and 460 in the 
case of a 4 x 4 x 4 supercell, were extracted for both spins, spanning the energy interval of 
(-14.0,9.1) eV with respect to the Fermi level (fully including the lowest valence bands). 
This energy interval is the outer window of WANNIER.90. 

At the Fermi surface, the relevant states are both dispersive conduction bands and 
hybridized cobalt states originating from the Co incomplete rf-shell. On the other hand, 
the closed copper (i-shell gives rise to narrow bands, lying deeper below the Fermi energy 
(/i). In order to simplify the host electronic structure, we want these Cu d-bands to 
be projected out during the disentanglement process and keep a simpler sp-like band. 
This can be achieved by both (1) the choice of the inner window position that freezes 
the bands with dispersive character around fi and (2) the choice of trial orbitals. As for 
the latter, we used one spherically symmetric function per atom, located in one of the 
two interstitials of the fee primitive cell. In addition, five cobalt centered orbitals with 
angular momentum I = 2 were used. The inner window choice is discussed in Section [3j 

§ For the supercell of 8 atoms (2x2x2), the Siesta calculation parameters are the following: the 
Brillouin zone sampling was 5x5x5, the discretization of the real-space mesh was taken from the 
mesh cutoff value 500 Ry and the self-consistency tolerance of the density matrix was 10~ 4 . 

For the supercell of 64 atoms (4x4x4), the Brillouin zone sampling was 3x3x3, keeping the 
same convergence criteria as for the 2 x 2 x 2-cell. 
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Figure 1. Scheme of the unit cell of the slab. Atomic basis functions are centered on 
Cu atoms (dots). Hence, the dots reproduce the slab geometry, where the size of the 
dots conveys information on the different x-coordinate. Crosses indicate the center for 
further Wannier.90 orbitals, the Wannier90 spherical interstitial trial orbitals, each 
located in the j and | of the bulk (111) translation vector, which points along the 
z-axis of the scheme. 



Finally, for the spread, Q, minimization a tolerance of 10 10 A 2 has been used. 
2.5. Simulations of Cu (111) 

The second system we studied is the Cu (111) surface. We used a 12-atom slab with 
a 1 x 1 surface unit cell. The k-point sampling was 14 x 14 x 1 and real space mesh 
cutoff of 800 Ry was used. Please, refer to [32] for more details. The Bloch functions 
corresponding to the lowest 120 bands were used to calculate the matrix M mn (k, b) 
using a k-point sampling 10 x 10 x 1 that largely suffices for convergence obtaining the 
MLWF. As in the previous case, we are interested in reproducing with the simplified 
scheme of MLWF the electronic structure about the Fermi energy, /i, hence we limit the 
inner energy window to the energy interval (—1.0,2.0) eV. The inner energy window 
must not contain more bands than the required number of Wannier functions, so the 
lower limit of the inner energy window is chosen somewhat above the narrow d-like 
bands and its upper edge of the inner window is limited by the higher-energy dispersive 
bands. 

Among various configurations of trial orbitals, we found that the one described in 
figure [1] converges to good-quality MLWF. The figure shows the centers of the Wannier 
functions that are located both at the atomic positions of the slab (dots) and in between 
(crosses) to enhance the electronic structure description. Inside the slab, one spherical 
interstitial is sufficient to give a good description of the electronic structure close to 
the Fermi energy, [i. On the surface, two functions are necessary; we put them at the 
positions corresponding to the interstitials of the bulk geometry. The trial orbital set has 
been found stable: only the outermost Wannier functions move some 0.21 A outwards 
during the spread minimization. 
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Table 1. Spread (f2), on-site energies (?) and occupation (n) for the cobalt Wannier 
functions of a substitutional Co atom in the 2 x 2 x 2-Cu cell. These split into two- 
and three-fold degenerate states m = {e g ,t2 g }. Majority (a =f) and minority (a 
spins are given. 
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Table 2. Spread (f2), on-site energies (?) and occupation (n.) for the cobalt Wannier 
functions of a substitutional Co atom in a 4 x 4 x 4-Cu cell. As in the previous case, the 
symmetry partially removes the degeneracy and now, two types of states are found, 
doubly and threefold degenerate: m = {e g ,t2 g } for the majority (a =t) and minority 
(a =1) spins. 
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0.67 



3. Results and discussion 

3.1. Cobalt impurity embedded in bulk copper 

In this section, we present and discuss MLWF for the impurity system. The essence of 
what we are doing is a tight-binding representation of a conduction band (supported 
by one interstitial Wannier function per atom) hybridizing with impurity Wannier 
functions. This picture will be extended in the next section, where we attempt to 
separate the Coulomb interaction in the impurity in the spirit of Anderson Hamiltonian. 

Figure |2] shows the ab-initio band structure and the disentangled one for the 2x2x2- 
cell. The color code denotes 100% overlap with a Co MLWF for red, and 0% for blue. 
The chosen inner energy window starts at —1.05 eV and extends up to +3.6 eV. This, 
in combination with trial orbital choice, efficiently disentangles the Cu sp-bands from 
the Cu d-bands, removes these last ones, and keeps the rest together with the full d- 
electron structure of the Co impurity. As we described above, the spin polarization of 
the system is sizable due to the presence of Co atoms. This is clearly seen in the present 
graph, where the bands split according to their spin, in figure [2] (a), the majority spin is 
represented, and correspond to Co <i-bands that coexist in the region of the Cu d-bands. 
The quality of these bands stemming from the MLWF calculation is considerably worse 
than the minority spin ones, figure |2] (b), because they lie outside the inner energy 
window. Nevertheless, the calculation contains information on the effect of the Cu d- 
bands on the majority spin Co bands as we will discuss in the analysis of the PDOS. 
Moreover, the calculations dealing with the magnetic structure and other information 
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close to the Fermi energy will be very accurate as can be seen in the excellent matching of 
the ab-initio bands and the MLWF ones near the Fermi level, figure [2j The calculation 
retains the Cu sp-bands that correspond to the colder- colored bands of figure [2] that 
strongly disperse. However, the hot-colored bands present flat features, revealing a 
small Co-Co intercell interaction, revealing that the dispersion largely comes from near- 
neighbor interactions, Co-Cu. For the 4x4x4 cell this interaction is even smaller, but 
the bands are not much flatter, what implies that for many Co-based properties our 
2x2x2 calculations will suffice. 

It can be seen that the MLWF disentanglement projects out the Cu <i-bands 
but retains a realistic description of Co d-bands, as revealed by looking into the spin 
polarization of the system. Indeed, by evaluating ([7]), we can compute the electronic 
population of the different MLWF. For the 4x4x4 cell, the majority spin MLWF yields 
4.59 electrons, while the minority spin gives a population of 2.65r;, if we add the spin 
polarization of the remaining Cu bands, we obtain 1.72 electrons, in good agreement 
with the 1.57 electrons of the full ab-initio calculation, discussed in section EU 

Valuable information is obtained by analyzing the PDOS in figures |3] and H] and the 
MLWF occupancies, on-site energies and real-space spreads as shown in tables [T] and EJ 

Due to the symmetry of the crystal field [22j [34], we find two types of MLWF, 
the e\ one at lower energy, which is twice degenerate according to the on-site energies 
(diagonal element of the MLWF Hamiltonian) , and the one of t 2g symmetry, three-fold 
degenerate. We find that the minority spin MLWF's, e^t^ are actually different in 
spread from the majority spin ones, ^ g -,t\ g . This is typical of unrestricted Hartree-Fock 
schemes, and it reflects the fact that both spins correspond to very different energy 
regions. The difference among MLWF's can be seen in the spread, Q in tables [U and 
[2] and in ([3]), that measures the extension of the MLWF. We see that e g MLWF's are 
more compact than the t 2g ones, and that the majority spin are more compact than the 
minority ones. It is interesting to note that the MLWF's are indeed very confined, in 
the present case, the more extended MLWF is basically zero beyond 2 A from its center. 

The on-site energy is exactly the first moment of the spectral function fl6]), and 
hence can be directly correlated to figures |3] and HJ The occupancy is the integration of 
the spectral function, over occupied states, from — oo to the Fermi energy, ([7]). Hence, 
these two quantities help us characterizing the spectral functions. As shown in figures [3] 
andSJ there are four sharp peaks that correspond to the above four types of MLWF. The 
presence of tails further from peak centers indicates that the MLWF's have an important 
hybridization with mainly the Cu host. These tails shift the on-site energy away from 
the peak energies. Most of the weight of the MLWF is, however, concentrated in a small 
energy window spanned by the peaks. Only in the case of the MLWF does the peak 
slightly lie outside the inner energy window used for the band disentanglement, and 
develops a tail in the Cu d-band region. 

Small differences show up when going from the 2x2x2 cell to the 4x4x4 one. 
The main peaks lie at the same energy, and only the first momenta or on-site energies 
are slightly different (table [2]) because the spectral function distribution is somewhat 
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Figure 2. (Color online) Interpolated band structure for bulk copper with a 
substitutional cobalt impurity in an eight-atom cell (2x2x2), (a) majority spin, 
(b) minority spin. Color denotes overlap of the eigenstate with cobalt states. For 
comparison, the ab-initio structure is plotted in grey. The zero of energy is the Fermi 
level. The inner energy window is the upper part of each diagram starting at the 
dashed line at —1 eV. 



narrower for the smaller cell. This is due to the smaller content of Co atoms in the large 
cell, that leads to smaller Co-Co hoppings. These comparison permits us to conclude 
that the 2x2x2 cell is indeed a good approximation to the dilute (paramagnetic) 
regime, in agreement with the experimental findings by Fan et al. [21] . 

The fact that the spectral functions are rather localized in energies permit us to 
conclude that they are very good descriptors of the actual Co electronic structure. This 
is further seen in figure |5l where the MLWF spectral function is compared with the 
atomic basis one (3d Cobalt PAOs). Indeed, the PAO spectral function is more spread 
in energies, and reveal more mixing with the Cu-band structure, and, consequently, less 
Co character. It is also interesting to notice that the main peaks coincide, showing that 
the MLWF contain a lot of physics of the actual electronic structure. 

From the spectral functions we see that the different electronic contributions of 
Co will have different physical properties. The t2 9 electronic structure is basically fully 
occupied and slightly contributes to magnetism. However the e g states have a large spin 
polarization, carrying the leading rhole in the magnetic properties of the Co impurities. 

3.2. Model Hamiltonian 

The results of the previous section show that the cobalt impurity physics can be 
qualitatively understood as an atomic <i-orbital weakly hybridized with the Cu substrate. 
The hybridization determines the positions and widths of atomic levels, which in turn 
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Figure 3. Cobalt Wannier-function projected density of states in the 4x4x4 cell 
for both spins. Zero energy coincides with Fermi level. 
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Figure 4. Cobalt Wannier-function projected density of states in the 4x4x4 cell 
for both spins. Zero energy coincides with Fermi level. 



fix the impurity occupancy. 

There is significant energy splitting regarding spin, which results in sizeable 
majority and minority spin populations. It is well understood that the origin of spin 
polarization lies in the intra-atomic Coulomb interaction [14] . However, the energy 
splitting for the e g levels is larger than for the t 2g ones, indicating that different t 2g and 
e g Coulomb matrix elements must be taken into account. 

The aforementioned physics can be captured by the multi-orbital Anderson 
Hamiltonian [H] 

(13) 

nk ma nka 
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Figure 5. Comparison of density of states projected on the cobalt e g and <2g Wannier 
functions and projected density of states onto I = 2 pseudo-atomic orbitals (PAO). 
The cell 4 x 4 x 4 is considered. Zero energy coincides with Fermi level. 



where the atomic part reads 



H a t ^ ^ £mJl"ma ^ ^ ^ U mm 'n m(T n m i CT , (1^0 



■i i i,i ii 



m,m =1,...5, a = — a. 

The substrate electrons (operators c^ kcr , c n k CT ) with Bloch energies e n k hybridize with 
impurity given by H at with on-site energies e m . Electrons on the impurity site are 
subject to Coulomb repulsion with matrix elements U mm i coupled to orbital occupancies 
n m a = d^dma. Interaction of equal spins is neglected by assuming the same value of 
direct and exchange integrals. The impurity and conduction band are connected by the 
hybridization term with matrix elements V n ^ m . 

Construction of the Hamiltonian ( ITS"]) from first principles using MLWF for a cobalt 
impurity in copper substrate is now discussed. Orthogonality of Wannier functions 
allows us to divide the Kohn-Sham Hamiltonian into the blocks 

(15) 

where H su bs is the Hamiltonian in the subspace of substrate Wannier functions, %i mp 
acts on the impurity Wannier functions (i. e. the cobalt e g and t2 9 functions) while the 
remaining terms V, are off-diagonal. In order to be consistent with the mean-field 
character of the Kohn-Sham DFT we identify (fl5|) with the mean-field solution of the 
Anderson Hamiltonian. 

Firstly, H at is discussed. We impose that the on-site energies e m obey restrictions 
of symmetry, i.e. have two values which will be denoted by e e and e t for e g and t 2g . In 
the correlation term, U mm i is a symmetric block matrix of the form 

(16) 




u ee 




Uet 


Utt 
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A usual mean-field factorization leads to the simplified Hamiltonian 

ma 

where the modified atomic energies read 

e eCT = e e + U ee n e „ + U et n w (18) 
tta = e* + Uunm + U et n e ^. (19) 

The right hand sides are identified with on-site energies of impurity Wannier functions, 
for we take H^j: F = Him P - The electron number n ecr means the total mean occupancy of 
e g symmetry orbitals with spin a and similarly for the t2 g states. 

To build the Anderson Hamiltonian, the bare energies e m and Coulomb integrals 
U mm > must be determined. Since the occupancies can be obtained from ([7]) we face a 
linear system of equations. 

This system is underdetermined. We solve it by fixing the difference 

A = e t - e e = 0.087e1/ (20) 

which is the crystal field splitting of atomic levels that is calculated as the t 2g — e g 
splitting of bulk copper d-orbitals . Taking the 4x4x4 case, we obtain: 

e t = -4.18eV, (21) 
U tt = 1.19eV, (22) 

(23) 

Eventhough we are not aware of any estimation of the U for cobalt in bulk copper, 
values reported in the literature for bulk Co are about 5 eV j?J , and for cobalt adsorbed 
on the (111) surface of gold [37] are 2.8 eV. We expect our values to be smaller than 
the ones obtainable by the approach of [7J. The reason for this is that our approach 
gives a rule to obtain an Anderson Hamiltonian from the MLWF Hamiltonian which 
proceeds from our DFT calculation and hence has the known problems of the local and 
semilocal exchange and correlation functionals in determining V (see for example [1]). 
Namely, our V will correspond to the Hund's rule exchange matrix element rather than 
the electrostatic one appearing in LDA+U. Nevertheless, our values albeit smaller, are 
comparable to other V values for Co. Indeed, Antonides and co-workers measured 1.2 
eV for the V of Co as an impurity [38] in excellent agreement with our calculation. We 
think that the main advantage of our approach is that it keeps the complete symmetry 
and multi-orbital character of the original DFT calculation. 

The construction of substrate and hybridization terms in (fTBl) is straightforward 
due to their one-particle form. Firstly, T-L su bs is diagonalized, yielding the conduction 

|| In good agreement with the crystal field splitting in Ag and Au estimated in [[35 J of about 0.15 eV, 
and in reference [36] of less than 0.1 eV. 



e e = -4.26eV, 
U ee = 1.63eV, 
Vet = 0.51eV. 
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Figure 6. Hybridization function, T, calculated for the cobalt impurity in the 2x2x2 
supercell of copper using a fine Brillouin zone sampling 40 x 40 x 40 and gaussian 
smearing 0.1 eV. The matrix elements corresponding to hybridization of e g states and 
t 2g states for both spins are given. Only diagonal contributions (see (|2~4"|) ) are shown 
because the non-diagonal hybridization functions are orders of magnitude smaller. Zero 
coincides with Fermi level. 



band energies e n k and Bloch states. In the next step, hoppings V and are transformed 
accordingly, leaving us with the matrix elements V n ^ m between substrate Bloch states 
and impurity Wannier functions. 

The essentially complete description of the impurity - substrate mixing is included 
in the hybridization function 



whose diagonal elements T mm give the spectral intensity of hybridization of a given 
impurity Wannier function m. The elements r mm (0) give (apart from a 27r factor) the 
inverse lifetime of an impurity electron. Off-diagonal terms provide substrate mediated 
intra-atomic hybridization intensity The function Y mm i{u) is precisely — - times the 
imaginary part of the retarded self-energy due to hybridization. 

So far we have tacitly omitted spin polarization inherent in the matrices % S ubst V 
and V', coming from a spin polarized Kohn-Sham DFT. The hybridization function 
calculated for majority and minority spin directions is given in figure |6j Off-diagonal 
£g ~ t2g matrix elements were omitted, for we found they are of the order of a few meV. 
The same holds for elements between different Wannier functions of the same symmetry. 
We see that all intensities have the same order of magnitude. The t 2g functions become 
much stronger for unoccupied substrate levels. The hybridization for minoritary spin 
is weaker than the majoritary-spin one. The most pronounced difference is between 
the e g t and e g 4- curves. The differences between hybridization functions for different 
spins partially reflect the polarization of the copper matrix by the magnetic moment of 
cobalt. 




(24) 
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Figure 7. Comparison of ab-initio band structure (grey) with the Wannier- 
interpolated (color) of a Cu (111) slab. Overlap of the eigenstate with the two 
outermost surface MLWF is indicated by color. The inner energy window spans the 
interval from -1.0 eV to 2.0 eV, where the matching between the interpolated and the 
ab-initio bands is very good. 

3.3. Surface studies: the case of Cu (111) 

Figure [7J shows the comparison of the ab-initio band structure for the Cu (111) and 
the MLWF interpolated ones. The disentangling scheme has permitted us to retain the 
electronic structure with sp character about the Fermi energy. As a consequence neither 
the d-bands nor the upper limit of the LL' gap (f for the surface Brillouin zone) are 
described by the MLWF. However, the electronic structure in the inner energy window 
is perfectly reproduced. 

In particular, we underline the excellent description of the Shockley surface state 
by the MLWF basis set. The band structure is basically indistinguishable from the 
ab-initio one which in turn is a very good description of the experimental one [32J. 
The minimization procedure leads to two types of differentiated MLWF sets, one set 
describing the bulk electronic structure with a spread of Q — 3.3 A 2 and a surface 
MLWF with n = 5.7 A 2 and its original center displaced by 0.21 A into the vacuum 
region. Hence, the MWLF try to follow the behavior of the PAO basis [32] where the 
energy minimization was improved by using two distinct basis sets, one for the bulk 
electronic structure and one for the surface with diffuse orbitals. Indeed, not only does 
the surface MLWF reproduce the surface state dispersion (minimum at E = —0.34 eV, 
and effective mass m* = 0.335 in electron masses, in good agreement with Eq = —0.42 
eV and m* = 0.37 of reference |32j) but it also shifts into the vacuum region in order to 
account for the surface spilling of charge. 

The interpolated bands can be analyzed in terms of the MLWF character they have. 
Figure [7J depicts in a color code the character of the interpolated bands. In this way we 
find that the Shockley surface state is basically purely described by the surface MLWF 
near F. As we move away from F, the surface state has some bulk MLWF weight, until 
it reaches the gap edge and it becomes a surface resonance with a large bulk MLWF 
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character. 

Outside the inner energy window, the bulk MLWF has more weight. Indeed, the 
lowest and highest energy bands are described by the bulk MLWF. This signals that 
in this slab geometry, the electronic structure has mainly a surface character and bulk 
MLWF are the least indicated to describe the slab electronic structure. 

The MLWF's correspond to s-like orbitals. A consequence of this is the 
disentanglement of the electronic structure as seen in figure [7J There, we see that 
the MLWF band structure spans from the bottom of the conduction band and crosses 
the d-band without mixing with it. Even outside the inner energy window, the MLWF- 
bands excellently reproduce the sp-bands. The description considerably worsens outside 
the interpolating inner energy window. 

4. Summary and conclusions 

We have implemented an interface to the scheme of [12J to obtain a finite set of 
maximally-localized Wannier functions (MLWF's) from calculations using the Siesta 
code [13] . 

The method yields an efficient band disentanglement in the case of Co impurities 
in bulk Cu. We have been able to retain just an s-like band to represent the Cu 
electronic structure, while keeping most of the correct description of the Co electronic 
structure. As a consequence, we have analyzed the magnetic properties in the limit of 
small concentration of Co impurities, and trace back most of the Co properties to the 
ones described by a single MLWF, doubly degenerated with e g symmetry. 

The reduction of the electronic structure to a few important elements is of great 
interest to obtain model Hamiltonians from full DFT calculations. We show that we can 
map the MLWF Hamiltonian into an Anderson model, having access to on-site energies, 
exact hybridization matrix elements and intra-atomic Coulomb matrix elements that 
reflect the correct symmetry of the problem. We show that by using this scheme, we 
obtain values of the intra-atomic Coulomb matrix that correspond to the semilocal 
exchange and correlation functional, and hence to the Hund's rule exchange, that are 
slightly smaller than the electrostatic ones obtained in LDA+U schemes jH [5j [6] . These 
results are encouraging for future studies of strongly correlated systems using MLWF. 

Even in the more stringent case of surface electronic structure, the description in 
terms of MLWF give very good results. We have applied the computational scheme to 
the Cu (111) surface and inside the chosen inner energy window the interpolation of the 
band structure is excellently reproducing the Shockley surface state with an accuracy 
comparable to the complete calculation of reference [32] , and permitting us to obtained 
the sp-band disentangled from the Cu d-bands. In this way, the electronic structure 
problem is limited to the sp-electronic structure about the Fermi energy. 

The surface case shows that not only is the MLWF method a mathematical trick 
to disentangle bands and reduce the problem to a smaller, tight-binding like, basis 
set. Indeed, one can analyze the obtained electronic structure in terms of the MLWF 
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and have interesting insight. We have projected the electronic bands in the two types 
MLWF of the surface problem: the surface and the bulk MLWF's. In this way, we have 
verified the mainly surface character of the Shockley state near the F point, as well as 
its subsequent mixing with the bulk MLWF as the surface state energy increases, until 
it becomes a surface resonance. 

In conclusion, we have implemented an interface to the scheme (12] to obtain a 
finite set of maximally-localized Wannier functions from calculations using the SIESTA 
code [13]. In this way, we can obtain DFT-based model Hamiltonians from ab-initio 
calculations with a very accurate description of the electronic structure inside a chosen 
energy window, that lends itself to the analysis and exploration of complex systems. 
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